__kernel void forward(int equationCount,
                      __global float* subDiagonal, __global float* mainDiagonal, __global const float* supDiagonal,
                      __global float* b)
{
    int blockIndex = get_global_id(0);
    
    int blockStartRow = blockIndex * ROWS_PER_BLOCK;
    int blockEndRow = blockStartRow + ROWS_PER_BLOCK - 1;

    float factor;

    for(int blockRowIndex = blockStartRow; blockRowIndex < blockEndRow; ++blockRowIndex)
    {
        factor = subDiagonal[blockRowIndex] / mainDiagonal[blockRowIndex];

        if(blockIndex > 0)
        {
            // Non-diagonal elements.
            subDiagonal[blockRowIndex] = -subDiagonal[blockRowIndex - 1] * factor;
        }
        else
        {
            subDiagonal[blockRowIndex] = 0.0;
        }

        mainDiagonal[blockRowIndex + 1] -= supDiagonal[blockRowIndex] * factor;
        b[blockRowIndex + 1] -= b[blockRowIndex] * factor;
    }
}

__kernel void backward(int equationCount,
                      __global float* subDiagonal, __global float* mainDiagonal, __global float* supDiagonal,
                      __global float* b)
{
    int blockIndex = get_global_id(0);

    int blockStartRow = blockIndex * ROWS_PER_BLOCK;
    int blockEndRow = blockStartRow + ROWS_PER_BLOCK - 1;
    blockStartRow = max(blockStartRow, 1);

    float factor;

    for(int blockRowIndex = blockEndRow - 1; blockRowIndex >= blockStartRow; --blockRowIndex)
    {
        factor = supDiagonal[blockRowIndex - 1] / mainDiagonal[blockRowIndex];

        if(blockIndex > 0)
        {
            if(blockRowIndex > blockStartRow)
            {
                subDiagonal[blockRowIndex - 2] -= subDiagonal[blockRowIndex - 1] * factor;
            }
            else
            {
                // If we edit element out of block, we take element from main diagonal.
                mainDiagonal[blockRowIndex - 1] -= subDiagonal[blockRowIndex - 1] * factor;
            }
        }

        // Non-diagonal elements.
        supDiagonal[blockRowIndex - 1] = -supDiagonal[blockRowIndex] * factor;

        b[blockRowIndex - 1] -= b[blockRowIndex] * factor;
    }
}

__kernel void auxiliaryMatrixCoeffs(int blockCount,
                                    __global const float* subDiagonal, __global const float* mainDiagonal, __global const float* supDiagonal,
                                    __global const float* b,
                                    __global float* auxiliaryMatrix)
{
    int blockIndex = get_global_id(0);

    int blockStartRow = blockIndex * ROWS_PER_BLOCK;
    int blockEndRow = blockStartRow + ROWS_PER_BLOCK - 1;

    if(blockIndex > 0)
    {
        auxiliaryMatrix[0 * blockCount + blockIndex - 1] = subDiagonal[blockEndRow - 1];
    }

    if(blockIndex < blockCount - 1)
    {
        auxiliaryMatrix[2 * blockCount + blockIndex] = supDiagonal[blockEndRow];
    }

    auxiliaryMatrix[1 * blockCount + blockIndex] = mainDiagonal[blockEndRow];
    auxiliaryMatrix[3 * blockCount + blockIndex] = b[blockEndRow];
}

__kernel void setX(__global const float* auxiliaryX,
                   __global float* x)
{
    int blockIndex = get_global_id(0);

    int blockStartRow = blockIndex * ROWS_PER_BLOCK;
    int blockEndRow = blockStartRow + ROWS_PER_BLOCK - 1;

    x[blockEndRow] = auxiliaryX[blockIndex];
}

__kernel void calculateX(int equationCount,
                         __global const float* subDiagonal, __global const float* mainDiagonal, __global const float* supDiagonal,
                         __global const float* b,
                         __global float* x)
{
    int blockIndex = get_global_id(0);

    int blockStartRow = blockIndex * ROWS_PER_BLOCK;
    int blockEndRow = blockStartRow + ROWS_PER_BLOCK - 1;

    // Use local variables and not arrays.
    float xEndRow = x[blockEndRow];
    float xStartRow;
    if(blockIndex > 0)
    {
        xStartRow = x[blockStartRow - 1];
    }

    float xRow;
    for(int blockRowIndex = blockStartRow; blockRowIndex < blockEndRow; ++blockRowIndex)
    {
        xRow = b[blockRowIndex];
        if(blockIndex > 0)
        {
            xRow -= subDiagonal[blockRowIndex - 1] * xStartRow;
        }
        xRow -= supDiagonal[blockRowIndex] * xEndRow;
        xRow /= mainDiagonal[blockRowIndex];

        x[blockRowIndex] = xRow;
    }
}
